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Entropic  lattice  Boltzmann  models  are  discrete- velocity  models  of  hydrodynamics 
that  possess  a  Lyapunov  function.  This  feature  makes  them  useful  as  nonlinear ly 
stable  numerical  methods  for  integrating  hydrodynamic  equations.  Over  the  last  few 
years,  such  models  have  been  successfully  developed  for  the  Navier-Stokes  equations 
in  two  and  three  dimensions,  and  have  been  proposed  as  a  new  category  of  subgrid 
model  of  tmrbulence.  In  the  present  work  we  develop  an  entropic  lattice  Boltzmann 
model  for  Burgers’s  equation  in  one  spatial  dimension.  In  addition  to  its  pedagogical 
value  as  a  simple  example  of  such  a  model,  our  result  is  actually  a  very  effective 
way  to  simulate  Biugers’s  equation  in  one  dimension.  At  moderate  to  high  values  of 
viscosity,  we  confirm  that  it  exhibits  no  trace  of  instability.  At  very  small  values  of 
viscosity,  however,  we  report  the  existence  of  oscillations  of  bounded  amphtude  in  the 
vicinity  of  the  shock,  where  gradient  scale  lengths  become  comparable  with  the  grid 
size.  As  the  viscosity  decreases,  the  amplitude  at  which  these  oscillations  saturate 
tends  to  increase.  This  indicates  that,  in  spite  of  their  nonlineair  stability,  entropic 
lattice  Boltzmann  models  may  become  inaccmate  when  the  ratio  of  gradient  scale 
length  to  grid  spacing  becomes  too  small.  Similar  inaccuracies  may  limit  the  utiUty 
of  the  entropic  lattice  Boltzmann  paradigm  as  a  subgrid  model  of  Navier-Stokes 
turbulence. 

Keywords:  lattice  Boltzmann;  Burgers’s  equation;  kinetic  theory 


1.  Introduction 

Lattice  Boltzmann  models  of  hydrodynamics  (Benzi  et  al.  1992;  Sued  2001)  evolve 
a  single-particle  distribution  function  in  discrete  time-steps  on  a  regular  spatiail  grid. 
The  velocity  space  is  discrete,  comprised  of  (possibly  hnear  combinations  of)  the 
lattice  vectors  themselves.  Though  this  constitutes  a  drastic  truncation  of  the  con¬ 
tinuum  Boltzmann  equation,  it  has  been  shown  that  the  derivative  hydrodynamic 
equations  may  emerge  unscathed  in  the  limit  of  small  Knudsen  number.  It  is  required 
only  that  the  kinetic  colhsion  operator  conserve  the  desired  hydrodynamic  quanti¬ 
ties,  such  as  mass,  momentum  and  energy,  and  that  the  lattice  possess  sufficient 
symmetry. 
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In  recent,  j'ears  it  has  been  shown  that  these  models  may  also  be  constructed  with  a 
Lyapunov  function — a.  discrete-time  and  discrete- velocity  version  of  Boltzmann’s  cel¬ 
ebrated  H  rlieoieui  as  is  appropriate  for  a  liydrf)flynaniic  ecpiatiou  with  di.ssijjatiou. 
In  this  case,  such  models  are  called  entropic  lattice  Boltzmann  models  (Boghosian  et 
al.  2001;  Karlin  et  al.  1998,  1999),  and  it  has  been  proposed  that  they  may  have  two 
other  very  desirable  features. 

(i)  First,  they  may  suppress  the  growth  of  numerical  instabilities,  and  thereby 
make  the  simulation  nonlineaidy  stable  (Ansumali  &  Karlin  2000,  2002;  Boghos¬ 
ian  et  al.  2001)  even  for  arbitrarily  small  viscosities.  This  argument  is  based  on 
the  premise  that  numerical  instabilities  evolve  in  ways  that  would  be  precluded 
by  the  existence  of  a  well-behaved  Lyapunov  function.  The  upper  hmit  to  the 
Reynolds  numbers  attainable  by  entropic  lattice  Boltzmann  models  is  therefore 
determined  by  loss  of  resolution  of  the  smallest  eddies,  rather  than  by  loss  of 
stabihty. 

(ii)  More  recently,  it  has  been  posited  that  entropic  lattice  Boltzmann  models  of 
the  Navier-Stokes  equations  constitute  a  natmal  subgrid  model  of  turbulence 
(Chen  et  al.  2003).  That  is,  even  if  the  grid  spacing  is  too  coarse  to  resolve  the 
smallest  eddies,  the  evolution  of  the  Lyaprmov  function  is  necessarily  mono¬ 
tonic,  and  it  is  argued  that  this  may  accurately  model  turbulent  dissipation  at 
subgrid  length-scales.  , 

Most  lattice  Boltzmann  models  employ  a  collision  operator  of  Bhatnagar-Gross- 
Krook  (BGK)  form,  with  collisional  relaxation  time  r  (Qian  et  al.  1992).  Since  trans¬ 
port  coefficients  such  as  viscosity  are  proportional  to  the  quantity  t  —  there  is 
incentive  to  make  r  as  small  as  possible  to  decrease  viscosity  and  thereby  increase 
Reynolds  number.  In  non-entropic  lattice  Boltzmann  models,  stability  is  guaranteed 
only  for  r  >  1,  severely  limiting  attainable  viscosities.  In  entropic  lattice  Boltz¬ 
mann  models,  by  contrast,  r  is  made  a  fimction  of  the  incoming  state  by  solving 
for  the  smallest  value  r*  <  1  that  does  not  increase  H.  The  value  then  used  is 
T  =  r*  -t-  k(1  —  T*),  where  0  <  /c  <  1  is  an  independent  variable.  It  has  been  shown 
that  the  viscosity  will  be  proportional  to  k,  and  may  therefore  be  made  arbitreuily 
small  by  the  appropriate  choice  of  k  (Boghosian  et  al.  2001;  Karhn  et  al.  1998,  1999). 
In  this  way,  the  entropic  lattice  Boltzmann  methodology  allows  for  arbitrarily  low 
viscosity  together  with  a  rigorous  discrete-time  H  theorem. 

In  earlier  papers,  we  constructed  entropic  lattice  Boltzmann  models  for  the  incom¬ 
pressible  Navier-Stokes  equations  that  are  Galilean  invariant  to  second  order  in  Mach 
nvunber,  and  we  showed  that  the  requirement  of  Galilean  invariance  makes  the  choice 
of  H  function  vmique  (Boghosian  et  al.  2003a,h).  More  specifically,  we  showed  that 
the  required  function  has  the  form  of  the  Burg  entropy  in  two  dimensions,  and  the 
Tsallis  entropy  in  higher  dimensions. 

In  this  work,  we  construct  an  entropic  lattice  Boltzmann  model  for  Burgers’s 
equation.  There  is  pedagogical  value  in  this,  since  it  provides  a  simple  illustration  of 
the  construction  of  such  a  model.  More  importantly,  however,  it  provides  a  simple 
model  with  which  to  test  the  above  hypotheses  about  entropic  lattice  Boltzmann 
models.  As  with  entropic  models  of  the  Navier-Stokes  equations,  we  may  make  the 
viscosity  arbitrarily  small,  giving  rise  to  shock  structures  with  large  gradients.  We 
may  thereby  test  whether  or  not  the  Lyapunov  function  truly  precludes  the  onset  of 
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numerical  instability,  and  we  may  investigate  the  utility  of  the  model  undei'  extreme 
conditions  when  shock  widths  become  comparable  with  the  grid  spacing. 

2.  Representation  of  the  distribution  function 

In  the  late  1980s,  it  was  shown  (Boghosian  &  Levermore  1987)  that  one-dimensional 
lattice-gas  models  with  a  conserved  mass,  a  maximum  occupation  number  per  site, 
and  a  wecik  spatial  asymmetry — that  is,  a  bias  for  movement  in  one  direction  over 
the  other — would  generically  give  rise  to  Burgers ’s  equation  hydrodynamics.  The 
magnitude  of  the  bias  should  scale  with  the  grid  spacing  in  the  continuum  limit.  Here 
we  recast  this  general  approach,  replacing  the  particulate  lattice-gas  model  with  a 
single-particle  distribution  function,  and  introducing  the  bias  in  the  H  function. 

We  consider  a  lattice  Boltzmann  model  on  a  one-dimensional  lattice,  with  a  two- 
component  distribution  function.  We  denote  by  N±  {x,  t)  the  number  of  particles  with 
velocity  ±1  at  site  x  and  time  t.  The  mass  at  position  x  and  time  t  is  then 

p{x,  t)  =  N+  {x,  t)  +  iV_  (x,  t) 

and  this  must  be  conserved  in  any  collision  process.  We  also  define  the  kinetic  degree 
of  freedom. 


u{x,  t)  =  iV+(ar,  t)  —  N_{x,  t), 

which  may  be  thought  of  as  a  non-conserved  momentum.  Together,  p  and  u  determine 
the  two  components  of  the  distribution  function: 

N±{x,t)  =  ^{p(x,t)  ±uix,t)).  (2.1) 

We  demand  that  the  components  of  the  distribution  function  be  positive  and  have 
an  upper  bound  of  unity: 

0  <  Al±  <  1. 

This,  in  turn,  implies  bounds  on  the  mass  and  momentum, 

0<  |(p±u)  ^  1, 

or 

0  ^  p  ^  2, 

max(— p,p  —  2)  ^  u  ^  min(-)-p,  2  —  p). 

3.  H  function  and  equilibrium 

Our  definition  of  the  H  function  must  be  such  as  to  allow  an  asymmetry  in  the 
collision  process.  We  suppose  that  the  H  function  is  of  trace  form,  and  weighted  by 
factors  a±  =  |(1  ±  a)  to  provide  the  asymmetry, 

H  =  a+h{N+)  -I-  a-h{N^), 


or 


l-)-a./p-|-u\ 

“H— ) 


-f- 


(3.1) 
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The  equilibrium  distribution  is  found  by  extreinizing  H  with  respect  to  u,  while 
keeping  p  fixed.  We  have 


(3.2) 


Note  that  if  a  =  0,  then  this  is  solved  by  u  =  0.  Moreover,  this  equation  is  invariant 
under  the  substitutions  u  — >  —u  and  a  — ^  —a.  It  follows  that  u  may  be  expanded 
as  a  perturbation  series  in  odd  powers  of  a.  To  do  this,  we  introduce  an  expansion 
parameter  e  by  making  the  formgd  replacement  a  — >  ea,  and  expanding  in  powers 
of  e\  at  the  end  of  the  calculation,  we  shall  reset  £  to  unity.  We  find 

-  ,  +  ■■■  , 


where 


,(eq,l)  _ 


-2a 


h'(\p) 

h"apy 


etc.  The  equilibrium  distribution  function  is  then  given  by 


N. 


(eq.o) 


-t-eN. 


(eq.l) 


Pe^N. 


(eq>2) 


+ 


where 


=  ip, 


;y^eq,2) 


h'{\p) 

h"apy 


(3.3) 


4.  Collision  operator 

To  perform  a  collision,  we  suppose  that  the  incoming  values  of  the  mass  and  velocity 
are  p  and  u.  We  compute  the  equilibrium  value  of  the  velocity  for  the  mass  p, 
and  then  we  implement  the  BGK  collision  operator  by  letting  the  post-collision  value 
of  u  be 

u' =  u  +  —  —  u).  (4.1) 

T 

We  choose 

r  =  r*  +  «(!  -  r*),  (4.2) 

where  k  6  [0, 1]  is  cui  input  parameter  that  shall  be  used  to  control  the  viscosity,  and 
T*  is  the  relaxation  time  that  keeps  H  fixed;  that  is 

h(^P,uP  -  u)^  =  H{p,u).  (4.3) 

Thus,  when  «  =  0,  the  collision  keeps  H  fixed,  and  H  is  a.  constant  of  the  motion. 
When  K  =  1,  the  collision  simply  replaces  the  incoming  state  with  the  equilibrimn. 
For  0<K<l,i/isa  Lyapimov  function,  always  decreasing  under  the  dynamics. 


Phil.  Vrans.  R.  Soc.  Land.  A  (2004) 


Entropic  lattice  Boltzmann  model  for  Burgers ’s  equation 


1695 


From  equations  (3.1)  and  (4.3)  we  get 


1  -I-  fV  /  f)  f  // 


+ 


1 


(“ 


I ) 


"j  + 


2r 


■{u 


III]) 


(4.4) 


This  must  be  solved  numerically  for  r*.  Such  solution  is  facilitated  by  the  certainty 
that  there  is  always  a  solution  t*  €  (|,  1)-  Equation  (4.1)  is  then  used  to  get  r, 
and  equation  (4.1)  then  gives  the  post- collision  value  of  u.  Finally,  this  is  used  in 
equation  (2.1)  to  get  the  post-collision  distribution  function. 

The  linearized  tollision  process  is  obtained  by  expanding  equation  (4.4),  treating 
a,  u  and  as  order  e,  and  simplifying  the  result  with  equation  (3.3).  We  find 
that  T*  =  \  +  0{e),  whence 

T  =  \  +  ^K-hO(e). 


Remarkably,  this  linearized  value  of  r  is  independent  of  the  incoming  state.  It  will 
be  useful  in  the  Chapman-Enskog  analysis  presented  in  the  next  section. 


5.  Chapman-Enskog  analysis 
The  lattice  BGK  equation  is 


N±{x±c,t  +  At)  -  N±(x,  t)  =  — -[iV±(x,  t)  —  f)], 

T 


and  this  may  be  rewritten  as 


|1  +  r 

which  has  the  formal  solution 


d 


exp(Af^±c 


^)-l  |iV±(a;,t)  =  7V(f‘’^(x,t), 


N±{ 


X,  t)  =  I 


1  -I-  r 


d 


exp  (  £  At—  ±  £c 


'dx 


-  1 


-1  OO 


j=0 


Here  we  have  introduced  parabolic  ordering  with  the  prescription  At  — >  e^At  and 
c  —>  ec.  The  distribution  fimction  is  thus  obtained  as  an  expansion, 

OO 


r=o 


where 


AtW  ^ 


=  N. 


(eq,0) 


(eq,2) 


T  TC 


dx 
dx 


—  r 


At 


^^£eq,0) 

dt 


8^ 
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We  immediately  obtain 

s-(0) 


A. 


il>- 


=  Ta 


Tcdp 


N. 


(2)  _ 


Tca 


h'{\p) 


d 


dx  [h"{\p)\ 


A 


ax-s 


fh'(ip)h"'(i/,) 

1  ^  _ 

T 

\^dp  Kc^d^p 

1  m^PW 

/  dx 

‘  2 

dt  2  9x2 

where  we  have  used  the  linearized  value  of  r  in  the  last  step. 

Finally,  we  demand  that  the  higher-order  corrections  to  the  distribution  function 
do  not  alter  the  definition  of  the  conserved  quantity  p  =  N+  +  1V_.  That  is,  we 
demand  that  -b  =  0  for  j  >  1.  We  see  that  this  is  automatic  for  j  =  1,  and 
that  for  j  =  2  we  get  the  hydrodynamic  equation 


dp  ca  ( h'(^p)h"'{^p)  }dp  _  d'^p 

r\  .  A.  I  r-.AA^1  r\  K  t  r\  ^ 


dt  At  I  [h"(ip)]2 


J  dx  2 At  dx'^ ' 


(5.1) 


6.  Recovery  of  Burgers’s  equation 


The  hydrodynamic  equation  obtained  in  the  last  section  will  be  of  the  form  of  Bmrg- 
ers’s  equation  if  we  demand  that  the  expression  in  brackets  in  equation  (5.1)  is  a 
linear  function  of  p.  To  this  end,  we  demand 

h'{z)h"'iz) 


or 


\h"{z)Y 

■  h'{z) 


1  =  —2az  —  h, 


2az  -t-  b, 


dz  \_h"{z) 

where  a  and  b  are  constants  to  be  chosen  later.  This  integrates  to  yield 

h'{z) 


h"{z) 


=  az^  bz  +  c, 


where  c  is  a  constant  of  integration.  If  we  define  g{z)  =  h'{z),  this  may  be  written 

^  1 

g{z)  az'^  +bz-]-  c 

It  is  simplest  to  choose  a  =  1  and  b—c=0,  so  that  this  integrates  to  yield 

g{z)  =  exp(-l/2) 

to  within  an  arbitrary  multiplicative  constant.  Other  choices  for  b  and  c  are  possible, 
and  may  yield  alternative  H  functions,  but  that  is  outside  the  scope  of  the  present 
work.  Finally,  we  integrate  one  more  time  to  obtain 

fi(z)  =  zexp(— l/z)  +  Ei(— 1/z),  (6.1) 
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Figure  1.  h{z)  versus  z. 


for  z  >  0,  where  ‘Ei’  is  the  exponential  integral  function.  This  monotonically  increas¬ 
ing  and  convex  function  is  plotted  in  figure  1. 

Using  the  form  of  h(z)  derived  in  equation  (6.1),  we  find  that  equation  (5.1)  reduces 
to 

dp  ca  dp  _  Kc^  d'^p 
dt  At^ dx  2At  dx'^ 

If  we  define  the  rescaled  dependent  variable 


ca 

w  =  --r-p, 
AV 


and  identify  the  viscosity, 


u  = 


KC? 


then  this  simplifies  to  Burgers’s  equation: 

dw  dvj  d^w 

Equation  (6.1)  for  h{z)  has  the  additional  virtue  of  allowing  for  an  analytic  solution 
for  the  equihbrium  distribution.  Inserting  it  into  equation  (3.2),  we  have 


(1  + 


Q!)exp(  — ^  ) 
\p  +  uj 


(1  —  a)  exp 


-2 


=  0, 


from  which  we  may  solve  exactly  for  the  equihbrium  value  of  u: 


2 

ln[(l  +  Q!)/(l  -  a)] 


(6.2) 


7.  Implementation 

We  implemented  the  above- described  entropic  lattice  Boltzmann  model  for  Burgers’s 
equation  on  grids  of  size  =  64  and  N  =  512,  with  initial  conditions 

p{x,  0)  =  po  +  pi  cos 
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We  took  a  =  —0.1,  po  =  0.8  and  pi  =  0.2  in  both  cases.  On  the  iV  =  64  grid  we  took 
K  =  0.0-5.  and  on  the  N  =  512  gi-id  we  took  k  =  0.01.  Both  of  these  were  intended 
as  extreme  ra.ses.  designed  ro  test  the  rohiistness  of  the  algorithm,  The  first  has  a 
very  coai'se  grid,  while  the  second  has  a  vei'y  low  viscosity  and  hence  a  very  nar  row 
shock.  In  both  cases,  the  ratio  of  shock  width  to  grid  spacing  is  small. 

A  regula  falsi  root  finder  was  used  to  solve  for  the  post-collision  velocity,  with  initial 
guess  -  u,  and  precautions  were  taken  so  that  the  iterates  never  left  the  feasible 
region  described  by  the  inequality  in  equation  (2.2).  Under  these  circumstfinces, 
convergence  was  very  robust.  It  is  likely  that  performance  could  be  improved  by  using 
the  regula  falsi  algorithm  to  get  in  the  vicinity  of  the  solution,  and  then  switching 
over  to  a  Newten-Raphson  solver,  but  we  did  not  explore  this  possibility. 

Figure  2  shows  the  N  =  64  simulation  for  time  increments  of  200,  from  t  =  0 
to  t  =  1800.  This  is  long  enough  to  see  the  formation,  steepening  and  decay  of  the 
shock.  By  t  =  1000,  it  is  clear  that  the  advsuicing  shock  is  leaving  small  oscillations 
behind  it.  These  oscillations  remain  bounded  in  amplitude,  and  eventually  decay 
to  zero  with  the  shock.  While  these  oscillations  are  not,  strictly  speaking,  a  Gibbs 
phenomenon — since  that  would  manifest  itself  both  ahead  of  and  behind  the  shock — 
they  are  similar  in  that  their  amplitude  remains  bounded,  and  goes  to  zero  as  the 
grid  is  refined  for  fixed  shock  width. 

Figure  3  shows  the  N  =  512  simulation  for  time  increments  of  1024,  from  t  =  0 
to  t  =  1800.  Once  again,  this  is  long  enough  to  see  the  formation,  steepening  and 
decay  of  the  shock.  By  t  =  4608,  it  is  clear  that  the  advancing  shock  is  leaving  small 
oscillations  behind  it.  This  oscillatory  tail  actually  grows  in  length,  while  remaining 
bounded  in  amplitude.  Once  again,  it  eventually  decays  to  zero  with  the  shock.  In 
both  cases,  if  the  rapid  oscillations  were  filtered  out,  the  remaining  smooth  compo¬ 
nent  of  the  solution  would  match  the  exact  solution  of  Burgers’s  equation  almost 
perfectly. 


8.  Conclusions 

We  have  derived  an  entropic  lattice  Boltzmemn  model  for  Bmrgers’s  equation,  and 
used  it  to  perform  a  fully  explicit,  imconditionally  stable  nmnerical  integration  of 
these  equations.  The  model  has  the  virtue  of  allowing  for  an  emalytic  solution  of  the 
equilibrium  distribution  function,  though  it  does  require  a  numerical  root-finding 
procedure  in  order  to  implement  the  collision  operator.  The  form  of  the  H  func¬ 
tion  that  is  required,  presented  in  equation  (6.1),  involves  the  exponential  integral 
function. 

We  have  also  used  this  model  to  simulate  Burgers ’s  equation  in  one  dimension  in 
two  extreme  situations:  the  first  with  an  unduly  coarse  grid,  and  the  second  with 
a  very  small  viscosity.  Both  situations  give  rise  to  small  ratios  of  shock  width  to 
grid  spacing.  In  both  situations,  we  noted  the  presence  of  oscillations  of  bounded 
amplitude  trailing  the  shock,  and  driven  by  the  strong  spatial  gradients  present 
there.  While  the  Lyapimov  function  did  indeed  prevent  the  unbounded  growth  of 
these  oscillations  into  a  full-fledged  numerical  instability,  their  presence  is  nonetheless 
troubling  if  one’s  goeil  is  to  use  entropic  collision  operators  as  subgrid  models  of 
turbulence.  Subgrid  turbulence  models  endeavour  to  do  ex£u:tly  what  we  have  done 
in  the  simulations  described  herein:  they  use  a  grid  spacing  that  is  known  to  be 
too  large  to  resolve  the  gradient  scale  lengths  present,  and  argue  that  the  eflfect  of 
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Figure  3.  p{x,  t)  versus  x  for  a  sequence  of  times  t  foi  N  =  512  and  k  =  0.01. 

tiorbulent  relaxation  at  finer  scales  is  captured  by  the  change  in  a  suitably  defined 
Lyapunov  function.  Thus,  entropic  models  used  for  this  purpose  may  likewise  be 
susceptible  to  oscillations  and  other  parasitic  dynamics  which,  though  bounded  in 
amplitude  and  scaling  to  zero  in  the  continuum  limit,  may  seriously  undermine  efforts 
to  use  the  model  for  the  accmrate  quantitative  study  of  turbulence. 

Having  said  the  above,  we  stop  short  of  suggesting  that  these  oscillations  preclude 
the  use  of  entropic  collision  operators  for  subgrid  studies  of  turbulence.  To  some 
extent,  the  tests  that  we  have  conducted  in  this  paper  are  unduly  stringent:  we 
have  noted  inaccuracies  in  the  evolution  of  the  model  for  one  particular  set  of  initial 
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conditions.  Most  siibgrid  models  endeavour  to  achieve  quantitative  faithfulness  only 
tp  quantities  that  can  be  averaged  robustly  over  ensembles  of  various  sorts.  For  this 
purpose,  iiiaccni  aeies  in  any  one  instantiation  may  be  iminiportant  a.s  long  as  they 
do  not  affect  the  ensemble  average.  Thus,  a  better  test  would  be  the  use  of  the 
model  described  herein  to  simulate  Burgers’s  turbulence  over  an  ensemble  of  initial 
conditions  on  ever  coarser  gi'ids.  We  leave  this  matter  to  future  work. 
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